#######################################
######Masculinity Threat Analysis######
#######################################
###############10/10/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science###
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis### 

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(purrr)
library(corrr)
library(psych)
library(corrplot)
library(forcats)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

MTData$P_EXP <- factor(MTData$P_EXP)
MTData$gender_binary <- factor(MTData$Q_CDGDRDESC_17, levels = c("1", "2"), labels = c("Men", "Women"))

dvs <- c("IRAQ_WAR_rescaled", "GAY_RIGHT_rescaled", "SUV_DESIRABILITY_rescaled", 
         "IMMIGRATION_rescaled", "ELECTRIC_CAR_rescaled", "TRANS_RIGHT_rescaled", "SUV_PAY_rescaled", 
         "HIRINGPOLICY_WOMEN_rescaled", "SYSTEM_JUSTIFICATION_rescaled",
         "MARIJUANA_rescaled", "CONSERVATISM_rescaled", "DOMINANCE_scale", "TRADITIONALISM_scale")

dv_labels <- c(
  IRAQ_WAR_rescaled            = "Iraq War",
  GAY_RIGHT_rescaled           = "Gay Rights",
  SUV_DESIRABILITY_rescaled    = "SUV Desirability",
  IMMIGRATION_rescaled         = "Immigration",
  ELECTRIC_CAR_rescaled        = "Electric Car",
  TRANS_RIGHT_rescaled         = "Trans Rights",
  SUV_PAY_rescaled             = "SUV Payment",
  HIRINGPOLICY_WOMEN_rescaled  = "Gender Affirmative Action",
  SYSTEM_JUSTIFICATION_rescaled = "System Justification",
  MARIJUANA_rescaled           = "Marijuana Legalization",
  CONSERVATISM_rescaled        = "Conservatism",
  DOMINANCE_scale              = "Social Dominance",
  TRADITIONALISM_scale         = "Traditionalism"
)

#######Correlations by Gender

plot_gender_correlations <- function(data, gender_label) {
  df <- data %>%
    filter(gender_binary == gender_label) %>%
    select(all_of(dvs)) %>%
    na.omit()
  
  # Correlation matrix with p-values
  cor_test <- psych::corr.test(df, use = "pairwise.complete.obs", method = "pearson", adjust = "none")
  
  # Apply labels
  colnames(cor_test$r) <- rownames(cor_test$r) <- dv_labels[colnames(cor_test$r)]
  colnames(cor_test$p) <- rownames(cor_test$p) <- dv_labels[colnames(cor_test$p)]
  
  # Plot
  corrplot::corrplot(cor_test$r,
                     p.mat = cor_test$p,
                     method = "color",
                     type = "upper",
                     tl.col = "black",
                     tl.srt = 45,
                     sig.level = 0.05,
                     insig = "blank",  # hide non-significant
                     mar = c(0, 0, 2, 0),
                     title = paste("Correlations among DVs:", gender_label))
}


pdf("correlations_women.pdf", width = 10, height = 10)
plot_gender_correlations(MTData, "Women")
dev.off()

pdf("correlations_men.pdf", width = 10, height = 10)
plot_gender_correlations(MTData, "Men")
dev.off()

#######Correlations by Experimental Group 

plot_condition_correlations <- function(data, condition_label) {
  df <- data %>%
    filter(P_EXP == condition_label) %>%
    select(all_of(dvs)) %>%
    na.omit()
  
  cor_test <- psych::corr.test(
    df,
    use = "pairwise.complete.obs",
    method = "pearson",
    adjust = "none"
  )
  
  # Apply cleaner labels
  colnames(cor_test$r) <- rownames(cor_test$r) <- dv_labels[colnames(cor_test$r)]
  colnames(cor_test$p) <- rownames(cor_test$p) <- dv_labels[colnames(cor_test$p)]
  
  corrplot::corrplot(
    cor_test$r,
    p.mat = cor_test$p,
    method = "color",
    type = "upper",
    tl.col = "black",
    tl.srt = 45,
    tl.cex = 0.7,
    sig.level = 0.05,
    insig = "blank",
    mar = c(0, 0, 2, 0),
    title = paste("Correlations among DVs:", condition_label)
  )
}

conditions <- c("no threat", "threat", "constrained", "PCKT")

# Display each plot on screen
for (cond in conditions) {
  cat("Plotting condition:", cond, "\n")
  plot_condition_correlations(MTData, cond)
}

# Save all plots to one PDF
pdf("correlation_plots_by_condition.pdf", width = 10, height = 10)
for (cond in conditions) {
  plot_condition_correlations(MTData, cond)
}
dev.off()


#######Correlation Between Trait Masculinity and Dependent Variables 

cor_results <- lapply(names(dv_labels), function(dv) {
  data.frame(
    DV = dv,
    correlation = cor(MTData[[dv]], MTData$THREAT_MALE_rescaled, use = "complete.obs")
  )
}) %>% bind_rows()

cor_results <- cor_results %>%
  mutate(DV_label = recode(DV, !!!dv_labels)) %>%
  mutate(DV_label = fct_reorder(DV_label, correlation))

ggplot(cor_results, aes(x = correlation, y = DV_label)) +
  geom_point(size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_minimal(base_size = 14) +
  labs(
    x = "Pearson Correlation with Trait Masculinity Threat",
    y = NULL
  )

ggsave(path = "Visualization", filename = "traitmasccorr.jpg", width = 15, height = 8, dpi=700)
dev.off()

#######Correlation Between Trait Femininity and Dependent Variables 

cor_results_fem <- lapply(names(dv_labels), function(dv) {
  data.frame(
    DV = dv,
    correlation = cor(MTData[[dv]], MTData$THREAT_FEMALE_rescaled, use = "complete.obs")
  )
}) %>% bind_rows()

cor_results_fem <- cor_results_fem %>%
  mutate(DV_label = recode(DV, !!!dv_labels)) %>%
  mutate(DV_label = fct_reorder(DV_label, correlation))

ggplot(cor_results_fem, aes(x = correlation, y = DV_label)) +
  geom_point(size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_minimal(base_size = 14) +
  labs(
    x = "Pearson Correlation with Trait Femininity Threat",
    y = NULL
  )

ggsave(path = "Visualization", filename = "traitfemcorr.jpg", width = 15, height = 8, dpi=700)
dev.off()